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We introduce a new geometric approach for the homogenization and 
inverse homogenization of the divergence form eUiptic operator with rough 
conductivity coefficients cr(x) in dimension two. We show that conduc- 
tivity coefficients are in one-to-one correspondence with divergence-free 
matrices and convex functions s{x) over the domain Q. Although homog- 
enization is a non-hnear and non-injective operator when apphed directly 
to conductivity coefficients, homogenization becomes a linear interpola- 
tion operator over triangulations of when re-expressed using convex 
functions, and is a volume averaging operator when re-expressed with 
divergence-free matrices. We explicitly give the transformations which 
map conductivity coefficients into divergence-free matrices and convex 
functions, as well as their respective inverses. Using optimal weighted De- 
launay triangulations for linearly interjjolating convex functions, we apply 
this geometric framework to obtain an optimally robust homogenization 
algorithm for arbitrary rough coefficients, extending the global optimal- 
ity of Delaunay triangulations with respect to a discrete Dirichlot energy 
to weighted Delaunay triangulations. Next, we consider inverse homog- 
enization, that is, the recovery of the microstructure from macroscopic 
information, a problem which is known to be both non-linear and severly 
ill-posed. We show how to decompose this reconstruction into a linear ill- 
posed problem and a well-posed non-linear problem. We apply this new 
geometric approach to Electrical Impedance Tomography (EIT) in dimen- 
sion two. It is known that the EIT problem admits at most one isotropic 
solution. If an isotropic solution exists, we show how to compute it from 
any conductivity having the same boundary Dirichlet-to-Neumann map. 
This is of practical importance since the EIT problem always admits a 
unique solution in the space of divergence-free matrices and is stable with 
respect to G-convergence in that space (this property fails for isotropic 
matrices). As such, we suggest that the space of convex functions is the 
natural space to use to parameterize solutions of the EIT problem. 
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1 Introduction 

In this paper, we introduce a new geometric framework of the homogenization 
and inverse homogenization of the divergence-form elliptic operator 

M^-div(CTVw) (1.1) 

where a is symmetric, uniformly elliptic and with entries S L°° . Owing to 
its physical interpretation, we refer to a as the conductivity. 

The classical theory of homogenization is based on abstract operator con- 
vergence and deals with the asymptotic limit of a sequence of operators of the 



form (1.1 1 parameterized by a small parameter e. We refer to G-convergence 
for symmetric operators, iJ-convergence for non-symmetric operators and F- 
convergence for variational problems [25J [321 HOI [SD [Ml [Zll IZS] . We also refer 
to [15] for the original formulation based on asymptotic analysis and [48 for a 
review. 

Instead of considering the homogenized limit of an e-family of operators of 
the form ( |1.1[ ) , we will construct in this paper a sequence of finite dimensional 
and low rank operators approximating ( |1.1[ ) with arbitrary bounded a{x). More 
precisely, since no small parameter e is introduced in the formulation, one has to 
understand homogenization in the context of finite dimensional approximation 
using a parameter h that represents a computational scale determined by the 
available computational power and the desired precision. We are motivated by 
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the fact that in most engineering problems, one has to deal with a given medium 
and not with an e-family of media. 

This observation gave rise to methods such as special finite element meth- 
ods, metric based upscaling and harmonic change of coordinates considered 
in |15l [TFl Uni [551 EZl EHl ISn] ■ This point of view recovers not only results from 
classical homogenization with periodic or ergodic coefficients but also allows for 
homogenization of a given medium with arbitrary rough coefficients. In partic- 
ular we need not make assumptions of ergodicity, scale separation and we do 
not need to introduce small parameter e. 

Our formalism, in not relying on small parameter e, is closely related to 
numerical homogenization which deals with coarse scale numerical approxima- 



tions of solutions of (2.11 below. Here we refer to the subspace projection 
formalism [65], the multiscale finite element method [45], the mixed multiscale 
finite element method [12], the heterogeneous multiscale method [36l|39], sparse 
chaos approximations [HI jBH] ; finite difference approximations based on viscos- 
ity solutions [57] , operator splitting methods [T^ and generalized finite element 
methods [75] . We refer to [37J[3H] for an numerical implementation of the idea 
of a global change of harmonic coordinates for porous media and reservoir mod- 
eling. 

Contributions. In this paper, we focus on the intrinsic geometric framework 
underlying homogenization. First we show that conductivities a can be put 
into one-to-one correspondence with, that is, can be parameterized by, sym- 
metric definite positive divergence free matrices Q, and by convex functions s 
as well (Section |2.2| . While the transformation which maps a into effective 
conductivities g'* per coarse edge element is a highly non-linear transformation 



(Section 2.11, we show that homogenization in the space of symmetric definite 
positive divergence free matrices Q acts as volume averaging, and hence is lin- 
ear, while homogenization in the space of convex functions s acts as a linear 
interpolation operator (Section [s]). Moreover, we show that homogenization 
as it is formulated here is self-consistent and satisfies a semi-group property 



(Section 3.3 1 



Hence, once formulated in the proper space, homogenization is a linear in- 
terpolation operator acting on convex functions. We apply this observation to 
construct optimally robust and stable algorithms for homogenizing divergence 
form equations with arbitrary rough coefficients by using weighted Delaunay tri- 
angulations for linearly interpolating convex functions (Section |4]). Figure [O] 
summarizes relationships between the different parameterizations for conduc- 
tivity we study. 

We use this new geometric framework for reducing the complexity of an in- 
verse homogenization problem (Section |5| . Inverse homogenization deals with 
the recovery of the physical conductivity a from coarse scale information, for in- 
stance, from the set of effective conductivities. This problem is ill-posed insofar 
as it has no unique solution, and the space of solutions is a highly nonlinear man- 
ifold. We use this new geometric framework to re-cast inverse homogenization 
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Figure 1.1: Relationships between parameterizations of conductivity. Straight 
and wavy hues represent linear and non-linear relationships, respectively. 
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into an optimization problem within a linear space. 

We apply this result to Electrical Impedance Tomography (EIT), the prob- 
lem of computing a from Dirichlet and Neumann data measured on the bound- 
ary of our domain. First, we provide a new method for solving EIT problems 
through parameterization via convex functions. Next we use this new geometric 
framework to obtain new theoretical results on the EIT problem (Section [6]). Al- 
though the EIT problem admits at most one isotropic solution, this isotropic so- 
lution may not exist if the boundary data have been measured on an anisotropic 
medium. We show that the EIT problem admits a unique solution in the space 
of divergence-free matrices. The uniqueness property has also been obtained 
in [5]. When conductivities are endowed with the topology of G-convergence 
the inverse conductivity problem is discontinuous when restricted to isotropic 
matrices [5l[53] and continuous when restricted to divergence-free matrices [5]. 

If an isotropic solution exists we show here how to compute it for any con- 
ductivity having the same boundary data. This is of practical importance since 
the medium to be recovered in a real application may not be isotropic and the 
associated EIT problem may not admit an isotropic solution but if such an 
isotropic solution exists it can be computed from the divergence-free solution 
by solving PDE (|6^. 

As such, we suggest that the space of divergence- free matrices parametrized 
by the space of convex functions is the natural space to look into for solutions 
of the EIT problem. 



2 Homogenization and parametrization of the 
conductivity space. 

To illustrate our new approach, we will consider, as a first example, the homog- 
enization of the Dirichlet problem 

f - div (aVu) = /, xen, 

[ u = 0, xedn. ^ ' ' 

S7 is a bounded convex subset of M'* with a boundary, and / G L°°(ri). The 
condition on / can be relaxed to / e L^(ri), but for the sake of simplicity we 
will restrict our presentation to / g L°°(ri). 



Let F : n ^ n denote the harmonic coordinates associated with (2.1 1. That 
is, F{x) = (Fi(a;), . . . , Fd{x)^ is a c?-dimensional vector field whose coordinates 
satisfy 

f div(CTVF,) = 0, xefl, 

S / (2-2) 

I Fi{x) = Xi, X G dil. 

In dimension d = 2 it is known that F is a homeomorphism from onto O and 
det(VF) > a.e. [6j [TOl [11]. For d > 3, F may be non-injective, even if a is 
smooth [inillUllS]. We will restrict our presentation to d = 2. 
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For a given symmetric matrix M , we denote by Aniax(-^^) and Aniin(-^^) its 
maximal and minimal eigenvalues. Define 



A,nax((V^^)^VF) 



(2.3) 



We will call condition (2.4 1 the following non-degeneracy condition on the 
anisotropy of {VF)'^VF 

< oo. (2.4) 



For d = 2, (2.4 1 is always satisfied if a is smooth [6] or even Holder continuous 

m- 

2.1 Homogenization as a non-linear operator 

Let ^Ih be a regular triangulation of having resolution h. Let Xh be the set of 
piecewise linear functions on 0,^ with Dirichlet boundary conditions. Let Afh be 
the set of interior nodes of fi/j. For each node i € Afh, denote ipi the piecewise 
linear nodal basis functions equal to 1 on the node i and on the other nodes. 
Let £h be the set of interior edges of O^, hence if e £ £h then e = (i, j) where i 
and j are distinct interior nodes and share the edge of two triangles of fi/j. 

2.1 Definition (Effective edge conductivities). Let q'^ be the mapping from £h 
onto M, such that for G £h 

I'^j -~ I ° F)fa{x)V{^, o F) dx. (2.5) 

Jn 

Observe that = q'^^, hence q'^ is only a function of undirected edges (i, j). We 
refer to q^^ as the effective conductivity of the edge (i, j). 

Let be the space of 2 x 2 uniformly elliptic, bounded and symmetric 
matrix fields on fi. Let T^h^^ be the operator mapping a onto q'' defined by 



(2.5 1. Let Qh be the image of T„h 



a — > Tqh^„[a] := q . 

Observe that T^h.o- is both non- linear and non-injective. 

Let j ~ i be the set of interior nodes j, distinct from i, that share an edge 
with i. 

2.2 Definition (Homogenized problem). Consider the vector {u^i)ieN'h of M-^'' 
such that for all i £ Afh, 

J2 - ) = / /(^)^^ ° n^) (2-7) 

We refer to this finite difference problem for {u^)i^j\/^ associated to q^ as the 
homogenized problem. 
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The identification of effective edge conductivities and the homogenized prob- 
lem is motivated by the following theorem: 



is unique. Moreover, let u be the solution of (2.1 ) and define 



2.3 Theorem. The homogenized problem (2.7 1 has a solution {u'^)ieA/h it 

(2.8) 



If condition (2.4) holds, then 



Uh\\Hi{n) < C/i||/||L~(n)- 



(2.9) 



2.4 Remark. We refer to 
theorem 12.31 



and 



for numerical results associated with 



2.5 Remark. The constant C depends on ||l/A,nin(a')||L=(o), || Aniax(f )||L°°(n)) 
Q, and fi. Replacing ||/||Loo(n) by ||/||L2(f2) in (2.9l adds a dependence of C on 
||(det(VF))-i||^^(^). 

2.6 Remark. Although the proof of the theorem shows a dependence of C on 



associated with condition (2.4 1, numerical results in dimension two indicate that 



C is mainly correlated with the contrast (minimal and maximal eigenvalues) of 
a. This is why we believe that there should be a way of proving (2.3 1 without 



condition (2.4). We refer to sections 2 and 3 of [B] for the detailed analysis of a 
similar condition (definition 2.1 of [B]). 



2.7 Remark. Problem (2.7) and Theorem 2.3 represent an generalization of 



method I of [16 to non-laminar media (see also [68 ). 

2.8 Remark. It is also proven in [68] (proof of Theorem 1.14) that if / G i°°(ri) 
then there exist constants C,a > such that u o F^^ E C^'"(ri) and 



||V(woF-i)||c«(j2) <C||/|Uoo(f,), 



(2.10) 



where constants C and a depend on fi, ||l/Amin(CT)||L~(o), || Aniax(o-)||L~(n), and 
/z. We also refer to [TB] (for quasi-laminar media) and [Bj for similar observations 
(on connections with quasi-regular and quasi-conformal mappings) 

2.9 Remark. Unlike a canonical finite element treatment, where we consider 
only approximation of the solution, here we are also considering approximation 
of the operator. This consideration is important, for example, in multi-grid 
solvers which rely on a set of operators which are self-consistent over a range of 
scales. 



Proof of Theorem \2.3\ The proof is similar to the proofs of Theorems 1.16 and 
1.23 of [68^; we also refer to [19^. For the sake of completeness we will recall its 
main lines. 
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Write Q the matrix (2.241. Replacing m by u o F in (2.1 1 we obtain after 
differentiation and change of variables that ii := uo F^^ satisfies 



(2.11) 



Similarly, multiplying (2.1 1 by test functions tpoF (with ip satisfying a Dirichlet 
boundary condition), integrating by parts and using the change variables y = 
F{x) we obtain that 



Jn Jn 



f 



det(Vi^) 



oF- 



(2.12) 



Observing that u satisfies the non-divergence form equation, we obtain, using 
Theorems 1.2.1 and 1.2.3 of [33], that if Q is uniformly elliptic and bounded, 
then u e W2'2(f]) with 



(2.13) 



The constant C depends on f2 and bounds on the minimal and maximal eigen- 
values of Q. We have used the fact that the Cordes-type condition on Q required 
by [Sn] simplifies for d = 2. 

Next, let Vfi be the linear space defined by (p o F for (p E Xh- Write Uh the 
finite element solution of (2.1 1 in Vh- Writing Uh as in (2.8 1 we obtain that the 
resulting finite-element linear system can be written as 

- q^^ - J2 9.' < = / /(^)^' ° ^(^) (2-14) 



for i e Afh- We use definition (2.5 1 for q'lj. Using the change of variables 
y = F{x) we obtain that g^v^- can be written 



4 



ipif" Q{x)'S/ (fij dx. 



Decomposing the constant function 1 over the basis ipj we obtain that 
- / {V^^fQ{x).v(}2vj)dx = Q, 



from which we deduce that 



(2.15) 



(2.16) 



(2.17) 



Combining (2.17) with (2.141 we obtain that the vector (uf)igjv'h satisfies equa- 
tion (|2J|. 
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Using the change of variables y = F{x) in 

{Viip^oF)fa{x).Vuhdx^ [ •froFf, (2.18) 
2 Jn 

we obtain that iih :— Uh o satisfies 

(V^.)-QV.. = /^^.^^o^-. (2.19) 

Hence uu is the finite element approximation of u. Using the notation <j[v\ := 
Jj^ Vw^trVw we obtain through the change of variables y — F{x) that a[v] = 
Q[vo F-^]. It follows that 

(t[u - Uh] = Q[u - Uh]- (2.20) 



Since Uh minimizes Q\u — v] over v G Xh we obtain equation (2.9) from the 
W2,2_j.ggyig^j.j^y ^ ( |2.13| ). □ 

The fact that g'', as a quadratic form on R-^'*, is positive definite can be 
obtained from the following proposition: 

2.10 Proposition. For all vectors {vi)itzj^/^ e 

E = / (^(^ ° i^))^aV(t; o F), (2.21) 

where v := Y.ieMu 

Proof. The proof follows from first observing that 

E"»9^>^ = I {^vf{y)Q{y){Vv){y)dy, (2.22) 

then applying the change of variables y — F{x). □ 

2.11 Remark. Despite the fact that positivity holds for any triangulation fi/j, 
as we shall examine in Section [4] we can take advantage of the freedom to 
choose Q.h to produce q^^ which give linear systems representing homogenized 
problems (2.7 1 having optimal conditioning properties. 



2.2 Parametrization of the conductivity space 

We now take advantage of special properties of a when transformed by its 
harmonic coordinates Q to parameterize the space of conductivities. 

2.12 Definition (Space of divergence-free matrices). We say that a matrix field 
M on is divergence- free if its columns are divergence-free vector fields. That 
is, M is divergence-free if for all v E and I E 

/ {Vvf M.l = (2.23) 
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2.13 Definition (Divergence-free conductivity). Given the conductivity a as- 
sociated to (1.1 1 and a domain fi, define Q to be the symmetric 2x2 matrix 
given by 



det(Vi^) 



(2.24) 



Again, F : ^ are the harmonic coordinates (2.2) associated to a. 



2.14 Proposition (Properties of Q). Q satisfies the following properties: 

• Q is positive-definite, symmetric and divergence-free. 

• det(Q) is uniformly bounded away from and oo. 

• Q is bounded and uniformly elliptic if and only if a satisfies the non- 
degeneracy condition (2.4 1. 



Prool Equations ( |2.1l| and ( |2.12| imply that for all u e iJ^ n H'^i^l) and all 

[ {V^pfQVu = - / ^y^Q.jd.djU. (2.25) 
Jn Jn 

Let Z e M'', choosing u such that Vw = I on the support of ip we obtain that for 
all leR'^ 

(VipfQ • / = (2.26) 



It follows by integration by parts that div{Q • Z) = in the weak sense and hence 
Q is divergence-free (its columns are divergence-free vector fields, this has also 
been obtained in [BS'). The second and third part of the Proposition can be 
obtained from 

det(Q) =det(croF-i), (2.27) 

and 

[ [ {VFfaWF (2.28) 
Jn Jn 

The last part of the Proposition can be obtained from the following inequalities 
(valid for d — 2). For x G a.e.. 



Amax('3) < A,nax(o')'( 

(0) > A„,i„(a). 



/ A.nax((Vf)^Vi^) 

A„,i„((VF)TVF) 
^ A^in((Vf)^VFy 

Amax((VF)^Vf^) 



(2.29) 

(2.30) 

□ 
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Write Tqjj the operator mapping a onto Q through equation (2.241. That 



IS, 



Tc 



M 

M 



{VFmVMVFm 
det(VFM) 



oF 



(2.31) 



M ' 



where Fm are the harmonic coordinates associated to M through equation (2.2 1 
(for CT = M) and A^div is the image of Ai under the operator Tq^o-- Observe 
(from Proposition 2.141 that Aldiv is a space of 2 x 2 of symmetric, positive and 



divergence-free matrix fields on il, with entries in L^(il) and with determinants 
uniformly bounded away from and infinity. 

Since for all M G A^div, Tq^AI] = M {Tq^^ is a non-linear projection 
onto A^div) it follows that Tq u is a non-injective operator from Ai onto A^div 
Now denote A^iso the space of 2 x 2 isotropic, uniformly elliptic, bounded and 
symmetric matrix fields on fi. Hence matrices in A^iso are of the form a{x)Id 
where Id is the d x d identity matrix and <j{x) is a scalar function. 

2.15 Theorem. The following statements hold in dimension d = 2: 

1. The operator Tq^^ is an injection from A^iso onto A^div 

2. 



To'AQ] = Vdet(Q)oG-i/d 



where G are the harmonic coordinates associated to 



Gi{x), j = 1, 2 is the unique solution of 



div 



ydet(Q) 



Gi{x) = Xi 



X eil, 
X e dn. 



(2.32) 

That is, 



(2.33) 



3. G = F ^ where G is the transformation defined by ( 2.33 1, and F are the 
harmonic coordinates associated to a :— Tq]^[Q] by (2.2 1. 



2.16 Remark. Observe that the non degeneracy condition ( 2.4 1 is not necessary 
for the validity of this theorem. 

2.17 Remark. Tq^^t is not surjective from A^iso onto A^div This can be proven 
by contradiction by assuming Q to be a non-isotropic constant matrix. Constant 
Q is trivially divergence-free, yet it follows that a — y^dct{Q) Id, F{x) — x and 
Q is isotropic, which is a contradiction. 

2.18 Remark. Tq ^ is not an injection from Ai onto A^dw • However it is 
known [GO. that for each a € A4 there exists a sequence in A^iso -ff-converging 
towards a. (Moreover, this sequence can be chosen to be of the form a(x, x/e), 
where a(x,y) is periodic in y.) Since A^iso is dense in Ai with respect to the 
topology induced by i?-convergence, and since Tg^o- is an injection from A^iso, 
the scope of applications associated with the existence of Tg ^ would not suffer 
from a restriction from Ai to A^iso- 
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Proof of Theorem 2.15 First observe that if a is scalar then we obtain from 

det(g) = {aoF-^f, 



equation (2.241 that 



and hence 



cr = A/det(Q) o F. 



(2.34) 



(2.35) 



Consider again equation (2.241. Let R be the 2x2, ^—rotation matrix in 
l^, that is, 

1 



R = 



(2.36) 



1 



det(A) 



RAR^. 



Observe that for a 2 x 2 matrix A, 

{A-r 

Write G := F'^ . Recall that 

VG = (VF)-ioi^-i. 



Applying (2.38) to (2.241 gives 



gVG = det(VG)((VG)~^)' aoG 



(2.37) 



(2.38) 



(2.39) 



Using ■y/det(Q) — a o G and applying equation (2.371 to ((VG) ^)'^ we obtain 
that 



Q 



-.VG^RVG R^. 



(2.40) 



yd^t(Q) 

Observing that in dimension two, for all functions v E H^, div(i?Vw) = 0, and 
we obtain from (2.401 that G satisfies equation (2.331. The boundary condition 
comes from the fact that G = F~^, where F is a diffeomorphism and F{x) — x 
on 90. Let us now show that equation (2.33 1 admits a unique solution. If G- is 
another solution of equation (2.331 then 

Q 



V(G. - G[) 



ydet(Q) 



V(G, - G:) = 



(2.41) 



Since Q is positive with entries and det(Q) is uniformly bounded away from 



zero and infinity it follows that 



is positive and its minimal eigenvalue is 



7dct(Q) 

bounded away from infinity almost everywhere in fi. It follows that V {Gi—G[) = 
almost everywhere in f2 and we conclude from the boundary condition on Gi 
and G\ that Gi = G[ almost everywhere in fi. □ 

2.19 Definition (The space of convex functions) . Consider the space of M^^'^(ri) 
convex functions on whose discriminants (determinant of the Hessian) are uni- 
formly bounded away from zero and infinity. Write S the quotient set on that 
space defined by the equivalence relation: s ~ s' if s — s' is an affine function. 
Let R be the rotation matrix (2.361. 
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2.20 Theorem (Scalar parameterization of conductivity in M^). For each Q € 
Aldiv there exists a unique s £ S such that 

Hess(s) = R^QR, (2.42) 

where Hess(s) is the Hessian of s. 

2.21 Remark. Since Q is positive-definite one concludes that Hess(s) is positive- 
definite, and thus, s{x) is convex. Furthermore, the principal curvature direc- 
tions of s{x) are the eigenvectors of Q, rotated by 7r/2. Note that this geometric 
characteristic will be crucial later when we approximate s{x) by piecewise-linear 
polynomials, which are not everywhere differentiable — but for which the notion 
of convexity is still well defined. 

Proof. In E^, the symmetry and divergence-free constraints on Q reduce the 
number of degrees of freedom of Q{x) to a single one. This remaining degree of 
freedom is s{x), our scalar convex parameterizing function. To construct s{x), 
observe that as a consequence of the Hodge decomposition, there exist functions 
G W^^^{n) such that 



o=(^ y^-i V 

where a, &, c are scalar functions. These choices ensure that the divergence-free 
condition is satisfied, namely that ax + hy = bx + Cy = 0. Another application 
of the Hodge decomposition gives the existence of s G W^'^iyi) such that Vs = 
{—k,h)^. This choice ensures that h — —hx — ky = —Sxy, the symmetry 
condition. The functions h and k are unique up to the addition of arbitrary 
constants, so s is unique up to the addition of affine functions of the type 
ax + Py + ^, where a, /?, 7 € K are arbitrary constants. □ 

We write T^.q the operator from A^div onto S mapping Q onto s. Observe 
that 

Ts,Q : Xdiv S 

is a bijection and 

T~q[s] = i?Hess(s)i?'^. (2.45) 

Refer to Figure |2.1| for a summary of the relationships between cr, Q and s. 
Figures [2T2] and [2T3| show an example conductivity in each of the three spaces. 



3 Discrete geometric homogenization 

We now apply the results of Section [2] to show that in our framework, homog- 
enization can be represented either as volume averaging, or as interpolation. 
Thus, unlike direct homogenization of a € M, homogenization in Mdiv or S is 
a linear operation. Moreover, in this framework, homogenization inherits the 
semi-group property enjoyed by volume averaging and interpolation, giving a 
self-consistency to homogenization in our setting. 
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Figure 2.1: The three parameterizations of conductivity, and the spaces to which 
each belongs. 




Figure 2.2: The left-hand image shows the original scalar conductivity a{x) = 
a(x) Id. In blue regions a=0.05, in red regions a=1.95, and in green regions, a 
= 1.0. The right-hand image gives ^det(Q) = cr o F~^, showing how harmonic 
coordinates distort (t(x). 
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Figure 2.3: Two views of the fine-scale function s{x) represented as a height 
field surface for the laminated conductivity of Figure |2.2| The left-hand view 
shows the fine-scale pattern in a{x), and the right-hand view highlights the 
coarse-scale anisotropy in the curvature. 



3.1 Homogenization by volume averaging 

The operator T^h u defined in (|2.6|) is a non- linear operator on A4. However, its 
restriction to A^divi which is a subset oi A4, is linear and equivalent to volume 
averaging as shown by Theorem |3.1| below. Using the notation of Section [271] 
we introduce the operator 



Tqh^ 



Mdiv — > Qh 



(3.1) 



where for Q G A^div and G £h, one has 

{T,.M)^, = - / (V<^O^QV^,. (3.2) 
Observe that T^k q is a volume averaging operator. 

3.1 Theorem (Homogenization by volume averaging). T^h q is a linear volume 
averaging operator on Mdiv Moreover: 



1. For Q & Mdiv one has 



2. For a eM 



(3.3) 
(3.4) 
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3. Writing xj the locations of the nodes of fi/i, for all I G 



(3.5) 



3.2 Remark. Equation (3.3 1 states that T^h q is the restriction of the operator 
to the space of divergence-free matrices A^div K follows from (3.4 1 that 
the homogenization operator T^h^^ is equal to the composition of the linear non- 
injective operator T^h q, which acts on divergence-free matrices, with the non- 
linear operator Tq u, which projects into the space of divergence-free matrices. 
Observe also that Tq^^ is injective as an operator from A^iso, the space of scalar 
conductivities, onto A^div 



3.3 Remark. Equation (3.5) is essentially stating that q is divergence free at 
a discrete level, see Section 2.1] for details. 

Proof. Using the change of coordinates y = F(x) we obtain that 



(3.6) 



which implies (3.4). One obtains equation (3.3) by observing that since Q is 
divergence-free its associated harmonic coordinates are just linear functions and 
T^^qIQ] = Q. Since Q is divergence-free, we have, for constant I e M?, 



{VLpifQ{x).ldx = 0. 



(3.7) 



Now, set Vh the set off all nodes in the triangulation Q,h and set Xj the location 
of node j £ Vh- The function z{x) SjeVh is the identity map on 

f2/j, so we can write I = ^ (^j(^i>^{l-Xj)ipj{x)j . Combining this with (3.7l 

gives (3.5). □ 



3.2 Homogenization by linear interpolation 

Write Tgh^g the linear interpolation operator over Uh- Hence for s £ S and 
:— Tsh g[s], we have, for x £ f2. 



s^{x) = y^j{xi)(pi{x), 



(3.8) 



where the sum in (3.8 1 is taken over all nodes of ilh and Xi is the location of 
node i. Write Sh the space of linear interpolations of elements of S on fi/j. 
Hence, 



s — > T^t^^sM ■■= ^ s{xi)(p^{x). 



(3.9) 
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Figure 3.1: Notation for computing q'^j from Si,Sj,Sk and si. 

For € £h write (5(j j)(a;) the uniform Lebesgue (Dirac) measure on the 
ed ge (as a subset of M^). Let R be the rotation matrix already introduced 



in (2.361. For G Sh observe that R}iess{s'^)R^ is a Dirac measure on edges 
oiQh- For (?, j) £ Sh define (T^h ,jh [s''])^^. as the curvature of in the direction 
orthogonal to the edge (z, j), hence 

T h h : Sh — > Qh 

, . 3.10 

with 

E (7^,N.45"]),/(»,j) = fiHess(s'')i?^. (3.11) 

Writ e Si — s{xi) where Si is the location of the node i. Referring to Fig- 
Tqfi sh[s''] defined as the trace of the i?-rotated Hessian of s'' on the 



ure 



3.1 



edgeTi.j) is 



{Tqh^^h [s'^])^. = - (cot e,jk + cot %i) Si 

- (cot ^jifc + cot Ojii) Sj (3.12) 

1 1 

Sfc + 7^, iSi, 



while diagonal elements (Tgh gh[s'^]) .^ are expressed as 
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where j ~ j is the set of vertices distinct from i and sharing an edge with vertex 
\eij\ is the length of edge \tijk\ is the area of the triangle with vertices 



(i, j, k), and 9ijk is the interior angle of triangle tijk at vertex j (see Figure 3.11 
Note that ( |3.12| ) is valid only for interior edges. Because of our choice to 
interpolate s{x) by piecewise linear functions, we have concentrated all of the 
curvature of s{x) on the edges of the mesh, and we need a complete hinge, 
an edge with two incident triangles, in order to approximate this curvature. 
Without values for s{x) outside of and hence exterior to the mesh, we do not 
have a complete hinge on boundary edges. This will become important where we 
apply our method to solve the inverse homogenization problem in EIT. However, 
for the homogenization problem, our homogeneous boundary conditions make 
irrelevant the values of on boundary edges. 

Tqh gh defined through (3.12) has several nice properties. For example, direct 



calculation shows that 



computed using (3.121 is divergence-free in the 



discrete sense given by (3.5 1 for any values s^. This fact allows us to parame- 



terize the space of edge conductivities g'' satisfying the discrete divergence-free 



condition (3.5 1 by linear interpolations of convex functions. 



3.4 Proposition (Discrete divergence- free parameterization of conductivity). 



Tqh gh defined using (3.121 has the following properties: 



1. Affine functions are exactly the nuUspace of T^h^^h ; in particular, 
T^fi j,h[s''] is divergence-free in the discrete sense of (3.5 1. 



2. The dimension of the range of T^U gh is equal to the number of edges in 
the triangulation, minus the discrete divergence-free constraints (3.5). 



3. Tqh gh defines a bijection from Sh onto and for s'' G Sh 



(3.13) 



Proof. These properties can be confirmed in both volume-averaged and inter- 
polation spaces: 



1. The first property can be verified directly from the hinge formula (3.12) 



For the Dirichlet problem in finite elements, the number of degrees of free- 
dom in a stiffness matrix which is not necessarily divergence-free equals 
the number of interior edges on the triangle mesh. The divergence-free con- 
straint imposes two constraints — one for each of the and y— directions — 
at each interior vertex such that the left term of ( 2.7 1 , namely ^ij {vi— 
Vj), is zero for affine functions. Thus, the divergence-free stiffness matrix 
has 



Ei - 2Vi 



(3.14) 



degrees of freedom, where Ej is the number of interior edges, and Vi is 
the number of interior vertices. 

The piecewise linear interpolation of s(x) has — 3 degrees of freedom, 
where there are V vertices in the mesh. The restriction of 3 degrees of 
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freedom corresponds to the arbitrary addition of affinc functions to s(x) 
bearing no change to Q. 

Our triangulation flh tessallates our domain fl, and so fl^ is a simply 
connected domain of trivial topology. For this topology, it can be shown 
that the number of edges E is 



E = 2V + Vi -3. 



(3.15) 



recalling that the number of boundary edges equals the number of bound- 
ary vertices, we have 

Ei -2Vi -3, (3.16) 

In fact s{x) and Q{x) represented on the same mesh have the same degrees 
of freedom when Q{x) is divergence-free. 



This property can be easily checked from the previous ones. 



□ 



3.5 Theorem. T^h^, a linear interpolation operator on S, has the following 
properties: 



1. For QeM div; 



2. For a eM, 



Tqh glQ] = Tqh ,.h O T^h ,, O Ts^q[Q]. 



(3.17) 



(3.18) 



3.6 Remark. It follows from equations (3.171 and (3.181 that homogeniza- 



tion is a linear interpolation operator acting on convex functions. Observe that 
jfi, s and T^ q are all linear operators. Hence, the non-linearity of the 
homogenization operator is confined to the non-linear projection operator Tg g. 
in (3.181 whereas if a is scalar its non-injectivity is confined to the linear in- 



terpolation operator Tgh ^. Equation (3.131 is understood in terms of measures 
on edges of and implies that the bijective operator mapping q'^ onto s'' is a 
restriction of the bijective operator mapping Q onto s to the spaces Qh and Sh- 

3.7 Remark. Provided that the Si interpolate a convex function s{x), the 
= (Tyh ,.h[s^fj^^ form a positive semi-definite stiffness matrix even if not all 

g^- are strictly positive. We discuss this further in the next section, where we 
show that even with this flexibility in the sign of the q'^j , it is always possible to 
triangulate a domain such that > 0. 

Proof. Define a coordinate system ^-77 such that edge ij is parallel to the r^-axis 



as illustrated in Figure 3.1 Using (2.151 to rewrite T^h qoTq^s[s] in this rotated 
coordinate system yields 



orjrj 



(3.19) 
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A change of variables confirms that integral (3.191 is invariant under rotation 
and translation. We abuse notation in that the second derivatives are un- 
derstood here in the sense of measures: we are about to interpolate s{x) by 
piecewise linear functions, which do not have pointwise second derivatives ev- 
erywhere. We are concerned with the values of s{x) interpolated at i,j,k, and 
Z, as these are associated to only the corresponding hat basis functions sharing 
support with those at i and j. The second derivatives of ip are non-zero only 
on edges, and due to the support of the gradients of the cp, contributions of the 
second derivatives at edges e^fc, ejk, en, and eji are also zero. Finally, the ds^nf 
and djirjifi are zero along ij, so the only contributions of s{x) to T^h q o Tq ^[s] 
defined through the integral are its second derivatives with respect to ^ along 
edge Cij . The contributions of four integrals remain, and by symmetry, we have 
only two integrals to compute. Noting that the singularities in the first and 
second derivatives are not coincident, from direct computation of the gradients 
of the basis functions and integration by parts we have 



dj^'fid^^ipidniP] = ^"-i^ + ' (3.20) 



d^ipid^^ipkd^ipj ^ - ^ (3.21) 

where \eij\ is the length of the edge with vertices (i,j), and \tijk\ is the area 
of the triangle with vertices {i,j,k). 9ijk is the interior angle of triangle ijk 



at vertex j (see Figure 3.11. The only contribution to these integrals is in the 



neighborhood of edge eij. Combining these results, we have that the elements 



of the stiffness matrix are given by formula (3.12). □ 



We refer to Figure [3T2| for a summary of the results of this subsection. 



3.3 Semi-group properties in geometric homogenization 

Consider, for example, three approximation scales < /ii < /12 < /13. We now 
show that homogenization from hi to /13 is identical to homogenization from hi 
to /12, then from /i2 to /13. We identify this as a semi-group property. 

Let ^Ihc be a coarse triangulation of fi, and ^Ihp be a finer, sub-triangulation 
of iliiF- Let ffjff be the piecewise linear nodal basis functions centered on 
the interior nodes of fi^^ and ^Ihp- Observe that for each interior node of the 
coarse triangulation i € Afhci ff can be written as a linear combination of (f-^, 
we write (jjtk the coefficients of that linear combination. Hence 

E '^^'^^fc- (3.22) 

Define T^hc ghp as the operator mapping the effective conductivities of the edges 
of fine triangulation onto the effective conductivities of the edges of the coarse 
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Physical 
conductivity space 
M (tensor) 
Miso (scalar) 



Divergence- free 
matrix space 



Convex functions 
space S 



(2.5), (2.71 



n, 




Hess(s) = BT^QR (2.421 



Q = RHess(s)R^ l |2.45| 



(3.91 



less{s^) = R'^q'^R \3.n\ r 



= RHess{s'^)R'^ i3.11 



effective edge 
conductivities 



Figure 3.2: Summary of discrete homogenization, showing the relationships be- 
tween the discrete spaces which approximate the spaces introduced in Section [2] 
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triangulation. Hence 

Tqhc „hp : Qhp Qhc 



(3-23) 



with, for G She 



(r,.c,, E Mmtf- (3-24) 

Let T^hp ^fep be the hnear interpolation operator mapping piecewise linear 
functions on Clhp onto piecewise linear functions on ^hc- Hence 

Tghc^ghp : Shp ^ She (3 25) 

s ^ > Tghc hp[s 



and as in (3.8 1, we have for x £ ft 

T,.a,s-As''ni^)^ E (3.26) 

3.8 Theorem (Semi-group properties in geometric homogenization) . The linear 
operators T^^c q''F E^nd T^h^ ^hp satisfy the following properties: 

1. Tghc ghp is the restriction of the interpolation operator T^hc ^ to piecewise 
linear functions on ilhp- That is, for s''^ g Shp 

T,.c,s-As''n = Ts^c,s[s''n- (3.27) 

2. For Q £ Mdw 

Tqi^c^giQ] = Tqhc^qhp o Tqhp glQ]. (3.28) 

3. For s e 5 

Tsi-cA^] = Tghc^s>-F ° T^'-F,sM- (3-29) 

4. For a e M 

'^q'^C,a[A = Tqhc^qhp O Tqhp^„[(7]. (3.30) 

5. For q^" e 

Tghc ,jhp [g''^] =T^hcr,s''c °T^'-c,s'^F ° ^q'^F^s'-F (3.31) 

6. For hi <h2 < /13 

^S^3 ,S^l — -^s'^3,S^2 ^ js'^l ■ (3.32) 

7. For hi < h2 < /13 

-^g''3,(j'»l ~ Tqh^ qh2 O Tqh2 qhi . (3.33) 
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Physical 
conductivity space 



Divergence-free 
matrix space 



Convex functions 
space 




Figure 3.3: Discrete geometric homogenization showing the sequence of scales 
referred to by the semi-group properties. 



3.9 Remark. As we will see below, if the triangulation flhj^ is not chosen 
properly s^'' — T^hp ^[s] may not be convex. In that situation T^h^ g in (3.271, 



when acting on s''^, has to be interpreted as a linear interpolation operator 
over f2/jp acting on continuous functions of fi. We will show in the next section 
how to choose the triangulation (resp. ^hc) to ensure the convexity of s^'' 
(resp. s'"^). 



3.10 Remark. The semi-group properties obtained in Theorem 3.8 are essential 
to the self-consistency of any homogenization theory. The fact that homogeniz- 
ing directly from scale hi to scale /13 is equivalent to homogenizing from scale 
hi to scale ft,2 then from /12 onto /13 is a property that is in general not satisfied 
by most numerical homogenization methods found in the literature when ap- 
plied to PDEs with arbitrary coefficients, such as non-periodic or non-ergodic 
conductivities. Figure [3?3| illustrates the sequence of scales referred to by these 
semi-group properties. 



4 Optimal meshes based on convex functions 

In this section, we use the convex function parameterization s S 5 to construct 
triangulations of which give matrices approximating the elliptic operator with 
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optimal conditioning. In particular, we show that we can triangulate a given 
set of vertices such that the off-diagonal elements of the stiffness matrix q^- are 
always non-positive. In turn, this minimizes the radii of the Gershgorin disks 
containing the eigenvalues of q^y The argument directly uses the geometry of 
s{x), constructing the triangulation from the convex hull of points projected 
up to s{x). We show that this procedure, a general case of the convex hull 
projection method for producing the Delaunay triangulation from a paraboloid, 
produces a weighted Delaunay triangulation. That is, we provide a geometric 
interpretation of the weighted Delaunay triangulation as well as an efficient 
method for producing optimal Q- adapted meshes. 

We remind the reader that throughout this paper, our triangulation Q.h is 
a tessellation of the compact and simply-connected domain fi, and, as such, is 
itself simply connected and of trivial topology. Also, since C M^, we shall 
identify the arguments of scalar functions as in s{x), x £ K^, or s{x, y),x,y £ K 
interchangeably without further comment. 



4.1 Construction of positive Dirichlet weights 



The constant C in (2.9 1 can be minimized by choosing the triangulation in a 



manner that ensures the positivity of the effective edge conductivities q^y The 



reason behind this observation lies in the fact that the discrete Dirichlet energy 



associated to the homogenized problem (2.7 1 is 



^E'?^'.("^-"^)' (4.1) 

where « ~ j are the edges of the triangulation, and Ui interpolate u{x) at vertices. 

We now show that for Q divergence-free, we can use a parameterization 
s{x) to build a triangulation such that > 0. g^ , identified here as Dirichlet 
weights are typically computed as elements of the stiffness matrix, where Q is 
known exactly. In this paper, we have introduced the parameterization s{x) 
for divergence-free conductivities, and if we interpolate s{x) by piecewise-linear 



functions, q^^ is given by the hinge formula ( |3.12 1 



In the special case where Q is the identity, it is well-known [70] that 

<1^J = ^(cot^.fe, +cot0,,,), (4.2) 

and in such case, all q^^ > when the vertices are connected by a Delaunay 
triangulation. Moreover, the Delaunay triangulation can be constructed geo- 
metrically. Starting with a set of vertices, the vertices are projected to the 
surface of any regular parabaloid 

p{x,y)^a{x^ + y^) , (4.3) 

where a > is constant. The convex hull of these points forms a triangulation 
over the surface of p(x,y), and the projection of this triangulation back to 
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the xy-p\ane is Dclaunay. See [M], for example. Our observation is that the 
correspondance 

Q = identity ^ s{x, v) ^ {x^ + V^) , (4.4) 

can be extended to all positive-definite and divergence-free Q. By constructing 
our triangulation as the projection of the convex hull of a set of points projected 
on to any convex s{x,y), we have the following: 

4.1 Theorem. Given a set of points V, there exists a triangulation of those 
points such that all g^- > 0. We refer to this triangulation as a Q-adapted 
triangulation. If there exists no edge for which g'^ = 0, this triangulation is 
unique. 

4.2 Remark. The points V D Vh, where Vh is the set of nodes in the re- 
sulting triangulation flh- That is, some points in V may be decimated by the 



triangulation. See also the remark following Proposition 4.5 



4.3 Remark, q^j > does not hold for arbitrary triangulations, as each trian- 
gulation, according to its connectivity, admits a different set of piecewise linear 
basis functions ipi. 

4.4 Remark. While s{x,y) may be convex, an arbitrary piecewise linear in- 
terpolation may not be. Figure [XT] illustrates two interpolations of s{x,y), one 
of which gives a q'lj > 0, and the other of which does not. Moreover, we note 
that as long as the function s(x, y) giving our interpolants Si is convex, the 
discrete Dirichlet operator is positive semi-definite, even if some individual ele- 



ments gf < 0. Figure 4.1 also illustrates how a Q-adapted triangulation can be 



non-unique: if four interpolants forming a hinge are co-planar, both diagonals 
give gf = 0. 



Proof of Theorem \4-l\ We proceed by constructing the triangulation as follows. 
Given V, we orthogonally project each 2D point onto the surface s{x) corre- 
sponding to Q. Take the convex hull of these points in 3D. Orient each convex 
hull normal so that it faces outward from the convex hull. Discard polyhedral 
faces of the convex hull with normals having positive z-components. Arbi- 
trarily triangulate polyhedra on the convex hull which are not already trian- 
gles. The resulting triangulation, once projected back orthogonally onto the 
plane, is the Q-adapted triangulation. Indeed, it is simple to show by di- 



rect calculation that hinge formula (3.12) is invariant under the transformation 



{si Si + aXi + byi + c}, where a, &, c G M are constants independent of i. This 
is consistent with the invariance of Q under the addition of affine functions to 
s{x). 

Now consider edge ij, referring to Figure [XT] Due to the invariance under 
affine addition, we can add the affine function which results in Si = Sj = Sk = 0. 
Thus, 

2\t,ji\ 
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Figure 4.1: Edge flips can replace non-convex edges, where Qab < 0, with 
convex edges without changing the interpolated values Si. For the given hinge, 
the diagonal giving a negative edge is on the left; a positive edge is on the right. 

\tiji\ is the unsigned triangle area, and so the sign of qij equals the sign of s;. 
That is, when s/ lies above the a;y-plane, qij > 0, showing that the hinge is 
convex if and only if qij > 0, and the hinge is flat if and only if = 0. All 
hinges on the convex hull of the interpolation of s(x, y) are convex or flat, so 
all qij > 0, as expected. Moreover, q^'^- = corresponds to a flat hinge, which in 
turn corresponds to an arbitrary triangulation of a polyhedron having four or 
more sides. This is the only manner in which the Q-adapted triangulation can 
be non-unique. □ 



4.2 Weighted Delaunay and Q-adapted triangulations 

There is a connection between s{x) and weighted Delaunay triangulations, the 
dual graphs of "power diagrams." Glickenstein [U] studies the discrete Dirichlct 



energy in context of weighted Delaunay triangulations. In the notation of (3.12 1 
and Figure |3.H Glickenstein shows that for weights Wi , the coefficients of the 
discrete Dirichlet energy are 

1i] = ^ (cot e^kj + cot Ouj) 

+ 12 (cot Oijk + cot 9iji) Wi 

(4.6) 

{cot 6 jik + cot 9 jii)wj 
1 



2 1 




1 








1 
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Comparison of this formula with (3.12) indicates that this is the discretization 
of 

9'; = - ^ (^2 - \Q^?j V(^„ (4.7) 

where I2 is the 2x2 identity matrix, and 



yy ^xy 

Wxy ^xx 

So, modulo addition of an arbitrary affine function, the interpolants 



(4.8) 



can be used to compute Delaunay weights from interpolants of s{x). 

Thus, we have demonstrated the following connection between weighted De- 
launay triangulations and Q-adapted triangulations: 

4.5 Proposition. Given a set of points V, the weighted Delaunay triangulation 
of those points having weights 

w, = xj + yj-2si (4.10) 

gives the same triangulation as that obtained by projecting the convex hull of 
points {xi,yi,Si) onto the a;?;-plane, where Si — s{xi,yi) are interpolants of the 
convex interpolation function s{x). 

4.6 Remark. Weighted Delaunay can be efficiently computed by current com- 
putational geometry tools, see for instance [2]. Thus, we use such a weighted De- 
launay algorithm instead of the convex hull construction to generate Q-adapted 
triangulations in our numerical tests below. 

4.7 Remark. In contrast to Delaunay meshes, weighted Delaunay triangula- 
tions do not necessarily contain all of the original points V. The "hidden" points 
correspond to values Si that lie inside the convex hull of the other interpolants 
of s{x,y). In our setting, as long as we construct Wi from Si interpolating a 
convex function s{x,y) (that is, weights representing a positive-definite Q), our 
weighted Delaunay triangulations do contain all the points in V. 

4.8 Remark. The triangulation is specific to Q, not to s{x,y). The addition 
of an affine function to s(x, y) does not alter the qij given by the hinge formula, 
a fact which can be confirmed by direct calculation. This is consistent with 
the observation that modifying the weights by the addition of an affine function 
{wi —> Wi + axi + byi + c}, a, 6, c € M are constants independent of i, does not 
change the weighted Delaunay triangulation. This can be seen by considering 
the dual graph determined by the points and their Delaunay weights, whereby 
adding an affine function to each of the weights only translates the dual graph 
in space, thereby leaving the triangulation unchanged. 
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4.9 Remark (Global energy minimum). The convex hull construction of a 
weighted Delaunay triangulation gives the global energy minimum result which 
is an extension of the result for the Delaunay triangulation. That is, the discrete 



Dirichlet energy (4.11 with q^j computed using hinge formula (3.121, where Si 
interpolate a convex s{x), gives the minimum energy for any given function 
Ui provided the 9^- are computed over the weighted Delaunay triangulation 



determined by weights (4.101 



To see this, consider the set of all triangulations of a fixed set of points. Each 
element of this set can be reached from every other element by performing a 
finite sequence of edge-flips. The local result is that if two triangulations differ 
only in a single flip of an edge, and the triangulation is weighted Delaunay after 
the flip, then the latter triangulation gives the smaller Dirichlet energy. 

A global result is not possible for general weighted Delaunay triangulations 
because the choice of weights can give points with non-positive dual areas, 
whereupon these points do not appear in the final triangulation. However, 
if the weights are computed from interpolation of a convex function, none of the 
points disappear, and the local result can be applied to arrive at the triangula- 
tion giving the global minimum of the Dirichlet norm. 

Similarly, if an arbitrary set of weights is used to construct interpolants Si 



using (4.10 1, taking the convex hull of these points removes exactly those points 
which give non-positive dual areas. See comments in |41] for further discussion 
of this global minimum result. 

4.3 Computing optimal meshes 

Using the connection that we established between s{x) and weighted Delaunay 
triangulations, we can design a numerical procedure to produce high quality 
Q-adapted meshes. Although limited to two-dimension, we extend the varia- 
tional approach to isotropic meshing presented in [9j to anisotropic meshes. In 
our case, we seek a mesh that produces a matrix associated to the homoge- 



nized problem (2.7 1 having a small condition number, while still providing good 
interpolations of the solution. 

The variational approach in |9| proceeds by moving points on a domain so as 
to improve triangulation quality. At each step, the strategy is to adjust points 
to minimize, for the current connectivity of the mesh, the cost function 

Ep= I b(x,y)-/(x,y)|, (4.11) 
Jn 

where p{x,y) — \{x'^ + y^) and p^{x,y) is the piecewise linear interpolation 
of p{x,y) at each of the points. That is, p^{x,y) inscribes p{x,y), and Ep 
represents the norm between the paraboloid and its piecewise linear inter- 
polation based on the current point positions and connectivity The variational 
approach proceeds by using the critical point of Ep to update point locations 
iteratively [3]. 
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Our extension consists of replacing the paraboloid p{x,y) with the conduc- 
tivity parameterization s(x,y). Computing the critical point of 



Es= / \six,y)-s\x,y)\ 
Jn 

with respect to point locations is found by solving 



(4.12) 



Hess(s)(x-,y*) 
Iless{s){x^,y^) 



^ s{xk ~ Xi,yk - yi) 



.ketj 



(4.13) 



Hess(s) is the Hessian of s(a;, y), Ki is the set 
is a triangle that belongs to Ki, and \tj \ is the 



for the new position {x* ,y*). 
of trinagles adjacent to point i, tj 
unsigned area of tj . Once the point positions have been updated in this fashion, 
we then recompute a new tessellation based on these points and the weights st 
through a weighted Delaunay algorithm as detailed in the previous section. 

4.10 Algorithm (Computing a Q-optimal mesh). Following [9 , our algorithm 
for producing triangulations that lead to well conditioned stiffness matrices for 



the homogenized problem (2.7 1 is as follows 



Read the interpolation function s{x) 
Generate initial vertex positions {xi,yi) inside 
Do 



Compute triangulation weights using (4.101 



Construct weighted Delaunay triangulation of the points 



Move points to their optimal positions using (4.131 
Until (convergence or max iteration) 



Figures [4~2] to [4~5] give the results of a numerical experiment illustrating the 
use of our algorithm for the case 



0.1 







10 



(4.14) 



Consistent with theory, the quality measures of interpolation and matrix con- 
dition number do not change at a greater rate than if an isotropic mesh is used 
with this conductivity. However the constants in the performance metrics of the 
anisotropic meshes are less than those for the isotropic meshes. 

As a word of explanation, for a fixed number of points on the boundary of 
the domain, anisotropic meshes tend to have fewer interior points than isotropic 
meshes. Since the experimental meshes are specified by their number of bound- 
ary points, this explains why the range of the total number of vertices is greater 
for the isotropic meshes than the anisotropic meshes. 
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Figure 4.2: Comparison of an isotropic and an anisotropic mesh. The figure on 
the left shows the lack of directional bias expected for a mesh suitable for the 
isotropic problem, while the figure on the right is suitable for the case where 
the conductivity is greater in the y-direction then in the x-direction. 




0.0001 I ' — ' — ' — ' — I 

10 100 1000 10000 100000 

number of vertices, N 

Figure 4.3: Interpolation quality of adapted meshes measured by the L2-norm 
error in a linear interpolation of s{x,y). Error diminishes as 0{N~^) in both 
cases, but is offset by a factor of about 4 in the adapted anisotropic meshes. 
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Figure 4.4: Interpolation quality of adapted mesiies measured by tiie /i^-semi- 
norm error in a linear interpolation of s{x,y). Error diminishes as 0{N~2) 
in both cases, but is offset by a factor of about 3 in the adapted anisotropic 
meshes. 
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Figure 4.5: Matrix conditioning quality of adapted meshes measured by the 
spectral condition number K2 of the stiffness matrix. The condition number 
grows as 0{N) in both cases, but is offset by a factor of about 5 in the adapted 
anisotropic meshes. 
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5 Relation with inverse homogenization 



Consider the following sequence of PDEs indexed by e. 

V / (5.1) 

u = 0, X G dil. 

Assuming that y a{y) is periodic (with in L°°(T''), where T'' is the d- 
dimensional unit torus) we know from classical homogenization theory [18 that 
Ue converges towards ug as e J, where uq is the solution of the following PDE 

f -div(CTeVuo) = /, xefl, 
Moreover ae is constant positive definite d x d matrix defined by 



a, := / aiy)ih + Vx{y)) (5.3) 

where the entries of the vector field x (xi; ■ • ■ i Xd) solutions of the cell 
problems 

-div(a(y)V(y, + x»(y))) = 0, y G 

(5.4) 



x^eH'{T^) / x^(.y) = o 

Jt'' 

Consider the following problem: 
Inverse homogenization problem: Given the effective matrix find a. 

This problem belongs to a class of problems in engineering called inverse 
homogenization, structural or shape optimization corresponding to the com- 
putation of the microstructure of a material from its effective or homogenized 
properties or the optimization of effective properties with respect to microstruc- 
tures belonging to an "admissible set" . These problems are known to be ill posed 
in the sense that they don't admit a solution but a minimizing sequence of de- 
signs. It is possible to characterize the limits of the sequences by following the 
theory of G-convergence 1^5\ as observed in |57l [62] . For non-symmetric ma- 
trices the notion of H-convergence has been introduced [STJ [55] . The modern 
theory for the optimal design of materials is the relaxation method through ho- 
mogenization l^ni HZl [5HI [55J [7S] . This theory has lead to numerical methods 
allowing for the design of nearly optimal micro-structures [U [TT] [5B] . We also 
refer to [30l [601 [HI] for the related theory of composite materials. 

In this paper we observe that at least for the conductivity problem in dimen- 
sion it is possible to transform the problem of looking for an optimal solution 
within a highly non-linear into the problem of looking for an optimal solution 



within a linear space, as illustrated by theorem 5.1 and Figure 5.1 for which 



efficient optimization algorithms could be developed. 



Define F{y) :~ y + x{v) where x is defined by (5.4 1 
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Figure 5.f : Illustration of theorem |5.1 



5.1 Theorem. Let Q be defined by (2.24), then 



1. Q is divergence-free, periodic and associated with a convex function s on 
[0,1] 2 through (|2.42|. 



2. 



3. If cr is isotropic then 



CTe = / Q{y) dy 



G - v/det(g) o G-i/d 
where G{y) ■.= y + x,X-^ (xi,X2) and 



(5.5) 
(5.6) 



— div 



Q 



^(y)v(2;. + Uy)) = 0, y e 

Vdet(Q) ) 



The problem of the computation of the microstructure of a material from 
macroscopic information is not limited to inverse homogenization. Indeed in 
many ill posed inverse problems one can choose a scale coarse enough at which 
the problem admits a unique solution. Hence these problems can be formulated 
as the composition of a well posed (eventually non-linear) problem with an in- 
verse homogenization problem. The approach proposed here can also used to 
transform these problems (looking for an optimal solution within a highly non- 
linear, non-convex space) into the problem of looking for an optimal solution 
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Figure 5.2: Relationships between spaces in inverse honiogenization. 



within a hnear space, as illustrated in Figure |5.2| for which efficient optimiza- 
tion algorithms can be used. As an example, we examine Electric Impedance 



Tomography (KIT) in subsection 6.1 



6 Electric Impedance Tomography 

We now apply our new approach to the inverse problem referred to as Elec- 
tric Impedance Tomography (EIT), which considers the electrical interpretation 



of (2.1 ). The goal is to determine electrical conductivity from boundary voltage 
and current measurements, whereupon a{x) is an image of the materials com- 
prising the domain. Boundary data is given as the Dirichlet-to-Neumann (DtN) 
map. Act : H^i^dil) H^^{d^), where this operator returns the electrical 
current pattern at the boundary for a given boundary potential. 
Act can be sampled by solving the Dirichlet problem 

I -div(CTVu) = 0, xe^, ^g^^ 
I u = g, X d dfl, 

and measuring the resulting Neumann data / = crf^,x e dfl. (In EIT, Neu- 
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mann data is interpreted as electric current.) Although boundary value prob- 



lem (6.11 is not identical to the basic problem (2.11, we can still appeal to 
the homogenization results in |68j provided we restrict g e W'^^p'^, in which 
case, we again have regular homogenization solutions u G W'^'^, that is, those 
obtained by applying conductivity Q{x) or s{x). If p > 2, a Sobolev embed- 



ding theorem gives u G C^'",a > 0, already seen in Section 2.1 although this 
restriction is not necessary for this section. 

The EIT problem was first identified in the mathematics literature in the 
seminal 1980 paper by Calderon |28j, although the technique had been known in 
geophysics since the 1930s. We refer to [2D] and references therein for simulated 
and experimental implementations of the method proposed by Calderon. From 
the work of Uhlmann, Sylvester, Kohn, Vogelius, Isakov and more recently, 
Alessandrini and Vessella, we know that complete knowledge of uniquely 
determines an isotropic a{x) e L°°i^), where ncR'',d>2 [71 liTl [5^ [75] . 

For a given diffeomorphism F from il onto fl, write 

It is known [33] (see also [THHH IMl [ZIl ) that for any diffeomorphism F : ft 
ri, F £ H^{il), the transformed conductivity d-{x) — F^,a{x) has the same 
DtN map as the original conductivity. If fj{x) is not isotropic, then a ^ a. 
However, as shown in [14 , this is the only manner in which a{x) € L°°(f2) can 
be non-unique. 

Let 5](ri) be set of uniformly elliptic and bounded conductivities on E M^, 
that is, 

^{ae L°°(r!;M2x2) \ a^a^,0< A„,i„(a) < X^M < oo}. (6.3) 

The main result of [TJ] is that A^- uniquely determines the equivalence class of 
conductivities 

e E(r!) I ai = 

(6.4) 

F : n n is an iJ -diffeomorphism and F \gfi— x}. 

It has also been shown that there exists at most 7 € i?o- such that 7 is isotropic [H] . 

Our contributions in this section are as follows: 
• Proposition |6.1| gives an alternate and very simple proof of the uniqueness 
of an isotropic a e E^-. This is in contrast to Lemma 3.1 of [14 , which 



appeals to quasi-conformal mappings. Moreover, Proposition (6.1 1 identifies 
isotropic 7 by explicit construction from an abribitrary M e E^. 
Proposition |6.2| shows that there exists equivalent classes E„ admitting no 
isotropic conductivities. 



Proposition 6.3 shows that a given a e E(ri) there exists a unique divergence- 



free matrix Q such that Aq — A„. It has been brought to our attention that 



proposition 6.3 has also been proven in 5 . Hence, although for a given DtN 



map there may not exist an isotropic a such that A = A^, there always exists 
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a unique divergence-free Q such that A — Aq . This is of practical importance 
since the medium to be recovered in a real application may not be isotropic and 
the associated EIT problem may not admit an isotropic solution. Although the 
inverse of the map a — > A^. is not continuous with respect to the topology of 
G-convergence when a is restricted to the set of isotropic matrices, it has also 
been shown in section 3 of that this inverse is continuous with respect to 
the topology of G-convergence when a is restricted to the set of divergence-free 
matrices. 

We suggest from the previous observations and from theorem |2.20| that 
the space of convex functions on f2 is a natural space in which to look for a 
parametrization of solutions of the EIT problem. In particular if an isotropic 



solution does exist, proposition (6.1 ) allows for its recovery through the resolu- 



tion the PDE (6.6 1 involving the hessian of that convex function. 



6.1 Proposition. Let 7 G i?cr such that 7 is isotropic. Then in dimension 
d= 2, 

1. 7 is unique. 

2. For any M £ E„ 



7= v/det(M)oG-i/d, 
the 

that is, G is the solution of 



(6.5) 

where G = (Gi, Go) are the harmonic coordinates associated to , , 

^ ' ^ ^dct(M)' 



div 



M 



v/det(M) 



VG,; = 0, 



X efi, 
X e dn. 



(6.6) 



3. G = F ^ where G is the transformation given by (6.6 1, and F is the 



diffeomorphism mapping 7 onto M through equation (6.2) 
Proof. The proof is identical to the proof of Theorem |2.15| 



□ 



6.2 Proposition. If cr is a non isotropic, symmetric, definite positive, constant 
2x2 matrix, then there exists no isotropic 7 G E^. 

Proof. Let us prove the proposi tion by contradiction. Assume that 7 exists. 



Then, it follows from Prop ositi on 6.1 that 7 is constant and equal to ^ydet{a)Id. 
Moreover, it follows from (6.6 1 that F^^{x) ~ x. Using 



(VF)^7V^ 
det(V^^) 

we obtain that a is isotropic which is a contradiction. 
6.3 Proposition. Let a £ E„. Then in dimension d=2, 



(6.7) 

□ 
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1. There exists a unique Q such that Q is a positive, symmetric divergence- 
free and Q — F^,a. Moreover, F are the harmonic coordinates associated 



to a given by (2.2) 



2. Q is bounded and uniformly elliptic if and only if the non-degeneracy 



condition (2.4 1 is satisfied for an arbitrary M ^ E„. 



3. Q is the unique positive, symmetric, and divergence-free matrix such that 

Proof. The existence of Q follows from Proposition |2.14[ Let us prove the 
uniqueness oi Q. If 

is divergence free, then for all I G and ip e C^{^) 

[ i^^fQ -1 = 0. (6.9) 
Jn 

Using the change of variables y — F{x) we obtain that 

[ {\/ipfQ-l^ f {\I{ipoF)faVF-l. (6.10) 

It follows that F are the harmonic coordinates associated to a which proves 
the uniqueness of Q. The second part of the proposition follows from Proposi- 
tion |2ll □ 

6.1 Numerical reconstructions with incomplete boundary 
measurements using geometric homogenization 

We close by examining two numerical methods for recovering conductivities 
from incomplete boundary data using the ideas of geometric homogenization. 
By incomplete we mean that potentials and currents are measured at only a 
finite number of points of the boundary of the domain. (For example, we have 



data at 8 points in Figure 6.2 for the medium shown in Figure 6.1 ) 

The first method is an iteration between the harmonic coordinates F(x) and 
the conductivity a{x). The second recovers s^{x) from incomplete boundary 
data, and from s^{x) we compute q'^{x), then Q. Both methods regularize the 
reconstruction in a natural way as to provide super-resolution of the conductivity 
in a sense we now make precise. 

The inverse conductivity problem with an imperfectly known boundary has 
also been considered in [SI] . We refer to [SD] and reference therein for an analysis 
of the reconstruction of realistic conductivities from noisy EIT data (using the 
D-bar method by studying its application to piecewise smooth conductivities). 
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Even with complete boundary data this inverse problem is ill-posed with 
respect to the resolution of a{x). The Lipschitz stability estimate in [7] states 



j^(l.W)(^)_^(2.A')(^)|| 



(6.11) 

where C{H^ {dfl), (dfl)) is the natural operator norm for the DtN map. 
(j^^'^\x) are scalar conductivities satisfying the ellipticity condition < A < 
(t(x) < almost everywhere in fl and belonging to a finite-dimensional space 
such that 

(6.12) 



1=1 



for known basis functions zl^\x). Thus, the inverse problem in this setting is 
to determine the real numbers crp'^^ from the given DtN map A^{j,N). 

The Lipschitz constant C{N) depends on and z'^^K As shown by con- 
struction in [71], when are characteristic functions of disjoint sets cov- 
ering C M'^, the bound 

C{N) > Aexp (^BN^^ (6.13) 

for absolute constants A,B>Ois sharp. That is, the amplification of error 
in the recovered conductivity with respect to boundary data error increases 
exponentially with N. 



From (6.131, we infer a resolution limit on the identification of a{x). Setting 
C our upper tolerance for the amplification of error in recovering a{x) with 
respect to boundary data error, and introducing resolution f = N~^/'^, which 
scales with length, we have 

- _ 

We refer to any features of <j{x) resolved at scales greater than this limit as 
stably-resolved and knowledge of features below this limit as super-resolved. 

6.1.1 Harmonic coordinate iteration. 

The first method provides super-resolution in two steps. First, we stably-resolve 
conductivity using a resistor-network interpretation. From this stable resolu- 
tion, we super-resolve conductivity by computing a function a{x) and its har- 
monic coordinates F{x) consistent with the stable resolution. 

To solve for the conductivity at a stable resolution, we consider a coarse 
triangulation of fl. Assigning a piecewise linear basis over the triangulation 
gives the edge-wise conductivies 



q-j - / {Vip^fQ{x)Vip,dx. ( |2l5 l 
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As we have already examined, when a{x) (hence Q{x)) is known, so too 
are the q^j. The discretized inverse problem is, given data at boundary vertices, 
determine an appropriate triangulation of the domain and the g'^- over the edges 
of the triangulation. We next specify our discrete model of conductivity in order 
to define what we mean by "boundary data." 

Let V/ be the set of interior vertices of a triangulation of fi, and let Vb be 
the boundary vertices, namely, the set of vertices on dil. Let the cardinality of 
Vb be Vb- Suppose vector u''^^ solves the matrix equation 

(6.15) 

ul ' \ I e Vb, 
where g^*^-* is given discrete Dirichlet data. Then we define 

/f^-E«^.K^'^-"f)' ^^^B (6.16) 

as the discrete Neumann data. The Vb linearly independent g^'^^ and their 
associated /'^'''^ together determine the matrix A^^, which we call the discrete 

Dirichlet-to-Neumann map. A^^ is linear, symmetric, and has the vector g = 

(1, 1, . . . , 1) as its nullspace. Hence, A^f has Vb{Vb — l)/2 degrees of freedom. 

In practise, the discrete DtN map is provided as problem data without a 
triangulation specified: only the boundary points where the Dirichlet and Neu- 
mann data are experimentally collected are given. We refer to this experimentally- 
determined discrete DtN map as A^^ . 

We are also aware that to make sense in the homogenization setting, the q^j 
must be discretely divergence-free. That is, we require that 

E qU^^^ - x?) = 0, * € V,, p = 1, 2, (6.17) 

where {xf'\x'f'^) is the ccy-location of vertex i. 

Set T^^ the set of triangulations having boundary vertices Vb specified by 
A^^, and {qij} the edge- values over T^^. The complete problem is 

minimise ||A\^ - A^«||*, 

r--A,i,} ' (6.18) 
subject to {Qj'^} discretely divergence-free. 

The norm \\-\\^ is a suitable matrix norm — as a form of regularization, we use a 
thresholded spectral norm which under-weights error in the modes associated to 
the smallest eigenvalues. We solve this non-convex constrained problem using 
Constrained Simulated Annealing (CSA). See [52], for example, for details on 
the CSA method. 
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EIT has already been cast in a similar form in [55], where edge-based data 
was solved for using a finite- volume treatment, interpreting edges of the graph 
which connects adjacent cells as electrical conductances. They determine the 
edge values using a direct calculation provided by the inverse theory for resistor 
networks [31] [35] . Although our work shares some similarities with this prior 
art, we do not assume that a connectivity is known a priori. 

An inversion algorithm for tomographic imaging of high contrast media 
based on a resistor network theory has also been introduced in [21]. The al- 
gorithm of [21] is based on the results of an asymptotic analysis of the forward 
problem showing that when the contrast of the conductivity is high, the current 
flow can be roughly approximated by that of a resistor network. Here our algo- 
rithm is based on geometric structures hidden in homogenization of divergence 
form elliptic equations with rough coefficients. 

Given an optimal triangulation T* and its associated stably-resolved {q^} 
representing conductivity, we now compute a fine-scale conductivity a-l'{x) con- 
sistent with our edge values, as well as its harmonic coordinates F{x). To help 
us super-resolve the conductivity, we also regularize a-^{x). 

Set a triangulation which is a refinement of triangulation T* from the 
solution to the stably-resolved problem. Let af{x) be constant on triangles of 
T-^ . Suppose coordinates F{x) are given, and solve 

minimise ||cr'^(x)||*, 

r f V. (6.19) 

subjectto -/(V(^,oF))^a/(a;)V(</7joF)=q,';., i,jer*. 

Here, ||-||* is some smoothness measure of {x). Following the success of regu- 
larization by total variation norms in other contexts, see jH |5T] for example (in 
particular we refer to (I2\ and references therein for convergence results the regu- 
larization of the inverse conductivity problem with discontinuous conductivities 
using total variation and the Mumford-Shah functional.), we choose 

||z(a;)||, = \\z{x)\\^v := [ |Vz(x)|. (6.20) 
Jn 

This norm makes particular sense for typical test cases, where the conductiv- 
ity takes on a small number of constant values. This "cartoon-like" scenario is 
common where a small blob of unusual material is included within a constant 



background material. The constraints in (6.19 1 are linear in the values of a-^ (x) 
on triangles of T-'^, and the norm is convex, so (6.191 is a convex optimiza- 
tion problem. In particular, it is possible to recast (6.191 as a linear program, 
see [51], for example. We then use the GNU Linear Programing Kit to do the 
optimization of this resulting linear program [3] . Note that we build our refined 
triangulation using Shewchuk's Triangle program |73j. 

The harmonic coordinates F{x) are not in general known. We set F(x) = x 



initially, and following the solution of (]6.19 1 , we compute 



-div(cr^VF =0, x£n, , , 

^ ^ ' ' (6.21 

F ^ X, X e dn, 
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Figure 6.1: A sample isotropic conductivity for testing reconstruction. The 
image on the left is cr, while the image on the right is \/SelQ = a o F^^. The 
dark blue background has conductivity 1.0, the red circle has conductivity 10.0 
and the yellow bar has conductivity 5.0. In this case, all of the features shrink 
in harmonic coordinates. 



using a' {x) from the previous step. We can now iterate, returning to solve (6.191 
with these new harmonic coordinates. 



Figures 6.1 and 6.2 show the results of a numerical experiment illustrating 
the method. In particular, the harmonic coordinate iteration resolves details of 
the true conductivity at scales below that of the coarse mesh used to resolve 
{<lij}- We observe numerically that this iteration can become unstable and fail 
to converge. However, before becoming unstable the algorithm indeed super- 
resolves the conductivity. We believe that this algorithm can be stabilized and 
plan to investigate its regularization in a future paper. 

6.1.2 Divergence-free parameterization recovery. 

Our second numerical method computes s{x) from boundary data in one step. In 
essence, we recover the divergence- free conductivity consistent with the bound- 
ary data, without concern for the fine-scale conductivity that gives rise to the 
coarse-scale anisotropy. 

We begin by tessellating by a fine-scale Dclaunay triangulation, and we 
parameterize conductivity by sj*, the piecewise linear interpolants of s{x) over 
vertices of the triangulation. From the s^, we can compute q^^ using the hinge 



formula in order to solve the discretized problem (6.15). This determines the 
discrete DtN map A^h in this setting. 

We shall also need a relationship between the and Q^^^j, an approximation 
of Q{x) constant on triangles. One choice is to presume that s{x) can be locally 
interpolated by a quadratic polynomial at the vertices of each triangle, and the 



opposite vertices of its three neighbours, see Figure 6.3 Taking second partial 



derivatives of this quadratic interpolant gives a linear relationship between Q'^^. 
and the six nearest . This quadratic interpolation presents a small difficulty, 
as triangles at the edge of Q. have at most two neighbours. Our solution is to 
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Figure 6.2: Output of the harmonic coordinate iteration. The figure on the 
top left is the coarse mesh produced by simulated annealing, the input to the 
harmonic coordinate iteration. Left to right, top to bottom, the remaining three 
images show the progression of the iteration at 1, 10, and 20 steps, showing its 
instability. The true conductivity is that of Figure |6.1[ 
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Figure 6.3: Stencil for approximating Q{x) on triangle ijk from the interpolants 
at nearby vertices. 

place ghost vertices outside the domain near each boundary edge, thus extending 
the domain of s{x) and adding points where must be determined. 

We solve for the sj* using optimization by an interior point method. Although 
the algorithm we choose is intended for convex optimization — the non-linear 
relationship between A^fi and the makes the resulting problem non-convex — 
we follow the practise in the EIT literature of relying on regularization to make 
the algorithm stable [23J|31j. We thus solve 

1 

minimise ;^ ^|| A^.g^'^) - + a||tr Q''||tv, ,^ 

^ k=i i^-^^) 

subject to Qij >0- 

We use IpOpt software package for the optimization [1]. 

The data are provided as K measured Dirichlet-Neumann pairs of data, 
{(gr('^\ /C^))}, and the Tikhonov parameter a is determined experimentally (a 
common method is the L-curve method). Again, the total variation norm is used 
to evaluate the smoothness of the conductivity. We could just as well regularize 
using det Q rather than trQ- Using the trace makes the problem more com- 
putationally tractable (the Jacobian is easier to compute), and our experience 
with such optimizations shows that regularizing with respect to the determi- 
nant does not improve our results. We compute the Jacobian of the objective's 
"quadratic" term using a primal-adjoint method, see [51], for example. 

We constrain q^j > on all edges, despite the possibility that our choice of 
triangulation may require that some edges should have negative values. Our 
reasons for this choice are practical: edge-flipping in this case de-stabilizes the 
interior point method. Moreover, numerical experiments using triangulations 
well-adapted to cr(x) do not give qualitatively better results. 



Figures pA] and [675] show reconstructions of the conductivities in Figures 6.1 



and]2.2[ respectively. We include the reconstruction of the conductivity in Fig- 



44 





Figure 6.4: Reconstruction of the isotropic conductivity in Figure 6.1 The 
left-hand figure shows trQ, while the right-hand figure shows -^/det Q. The 
reconstruction blurs the original a, similar to other methods in the literature, 
but does not underestimate the dynamic range of the large rectangle. 



ure |6.1| only to show that our parameterization can resolve this test case, a 
typical one in the EIT literature. For such tests recovering "cartoon blobs," 
our method does not compete with existing methods such as variational ap- 
proaches [23], or those based on quasi-conformal mappings [ISl ^9j. Our re- 
covery of the conductivity in Figure |2.2| however, achieves a reconstruction, 
to our knowledge, which has not previously been realized. The pitch of the 
laminations in this test case are below a reasonable limit of stable resolution. 
Hencee do not aim to recover the laminations themselves, but we do recover 
their up-scaled representation. The anisotropy of this up-scaled representation 
is apparent in Figure |6.5[ Admitting the possibility of recovering anisotropic, 
though divergence-free, conductivities by parameterizing conductivity by s{x) 
gives a sensical recovered conductivity. Owing to the stable resolution limit, 
parameterizing (j{x) directly, by a usual parameterization — such as the linear 
combination in ( 6.12[ ), choosing Zi{x) as characteristic functions — is not success- 
ful. 
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